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Abstract. Biochemical models that exhibit bistability are of interest 
to biologists and mathematicians alike. Chemical reaction network the- 
ory can provide sufficient conditions for the existence of bistability, and 
on the other hand can rule out the possibility of multiple steady states. 
Understanding small networks is important because the existence of mul- 
tiple steady states in a subnetwork of a biochemical model can sometimes 
be lifted to establish multistationarity in the larger network. This paper 
establishes the smallest reversible, mass-preserving network that admits 
bistability and determines the semi-algebraic set of parameters for which 
more than one steady state exists. 
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1 Introduction 

Bistable biochemical models are often presented as the possible underpinnings 
of chemical switches [2117] . Systematic study of mass-action kinetics models- 
which a priori may or may not admit multiple steady states-constitutes chemical 
reaction network theory (CRNT), a subject pioneered by Horn, Jackson, and 
Feinberg [13I16| . Certain classes of networks, such as those of deficiency zero, 
do not exhibit multistationarity or other strange behaviors. A generalization 
of deficiency-zero systems is the class of toric dynamical systems which have a 
unique steady state [S]. See also the recent work of Craciun and Feinberg for 
additional conditions that rule out multistationarity [617] . 

On the other hand, there are conditions that are sufficient for establishing 
whether a network supports multiple steady states. The CRNT Toolbox de- 
veloped by Feinberg and improved by Ellison implements the Deficiency One 
and Advanced Deficiency Algorithms [9112] ; this software is available online [10] . 
For a large class of systems, the CRNT Toolbox either provides a witness for 
multistationarity or concludes that it is impossible. For systems for which the 
CRNT Toolbox is inconclusive, see the approach of Conradi et al. [I]. Related 
work includes an algebraic approach that determines the full set of parameters 
for which a system is multistationary; a necessary and sufficient condition for 
multistationarity is the existence of a non-trivial sign vector in the intersection 
of two subsets of Euclidean space [5] . 



To model biological processes, one typically reverse-engineers a system of 
non-linear differential equations that exhibits specific dynamical behavior, such 
as bistability or oscillations, observed in experiments. For example, Segel pro- 
poses a small immune network consisting only two cell types, which has three 
stable steady states, corresponding to "normal," "vaccinated," and "diseased" 
states |f 9j . Similarly, the Brusselator is a mass-action kinetics network with a 
stable limit cycle [TTj . 

This paper focuses on the smallest mass-action kinetics networks that admit 
multiple steady states. Section 2 provides an introduction to chemical reaction 
network theory. A special network called the Square is shown in Section 3 to be 
a smallest reversible multistationary chemical reaction network. Sections 4 and 5 
determine precisely which parameters of the Square give rise to multiple steady 
states. 

2 Chemical Reaction Network Theory 

We now give an introduction to chemical reaction network theory. Before giving 
precise definitions, we present an intuitive example that illustrates how a chemi- 
cal reaction network gives rise to a dynamical system. An example of a chemical 
reaction, as it usually appears in the literature, is the following: 

A + B ^—^ 3A + C 

In this reaction, one unit of chemical species A and one of B react (at reaction 
rate k) to form three units of A and one of C. The concentrations ca, Cb, and cc 
will change in time as the reaction occurs. Under the assumption of mass-action 
kinetics, species A and B react at a rate proportional to the product of their 
concentrations, where the proportionality constant is the rate constant k. Noting 
that the reaction yields a net change of two in the amount of A, we obtain the 
first differential equation in the following system: 

d 

—ca = 2kcaCb , 
at 

d 

— C B = - KC A C B , 

at 
d 

— C C = KCACB ■ 

dt 

The other two equations arise similarly. Next we include the reverse reaction 
and switch from additive to multiplicative notation to highlight the monomials 
that appear in our differential equations; the chemical reaction networks in this 
paper will appear with the following notation: 



K 

cacb * c\c c 



This network defines differential equations that are each a sum of the mono- 
mial contribution from the reactant of each chemical reaction in the network: 




The recipe for obtaining these differential equations from a chemical reaction 
network easily generalizes from this example. However, in order to display the 
linearity hidden in these non-linear equations, the equations will appear in a 
different but equivalent form in (fT]) below. 

We now establish the notation for this paper, following 5J. A chemical reac- 
tion network is a finite directed graph whose vertices are labeled by monomials 
and whose edges are labeled by parameters. Specifically, the digraph is denoted 
G = (V, E), with vertex set V = {1, 2, . . . , n} and edge set E C {(£, j) e V X V : 
i j}. The vertex i of G represents the ith chemical complex and is labeled by 
the monomial 

r Vi — r Vil r Vi2 , . , y is 

This yields Y — (ytj), an n x s-matrix of non-negative integers. The unknowns 
c\ , C2, . . . , c s represent the concentrations of the s species in the network, and we 
regard them as functions Ci{t) of time t. The monomial labels form the entries 
in the following row vector: 



A network is said to be mass-preserving if all monomials c Vi have the same 
degree. Each directed edge (i, j) £ E is labeled by a positive parameter Kij which 
represents the rate constant in the reaction from the i-th chemical complex to 
the j-th chemical complex. A network is reversible if the graph G is undirected, 
in which case each undirected edge has two labels Kij and Kji. Let A K denote the 
negative of the Laplacian of the digraph G. In other words A K is the n x n-matrix 
whose off-diagonal entries are the Kij and whose row sums are zero. Mass-action 
kinetics specified by the digraph G is the dynamical system defined by 



By decomposing the mass-action equations in this way, we see that they are linear 
in the Kij by way of the matrix A K . A steady state (or equilibrium) is a positive 
concentration vector c S M^q at which the equations (TTJ vanish. These equations 
remain in the (stoichiometric) subspace S spanned by the vectors yi — yj (where 
(i, j) is an edge of G). In the earlier example, yi — yj = (—2, —1, 1), meaning 
that whenever a reaction occurs, two units of A and one of B are lost, while one 
unit of C is formed (or vice- versa). Therefore, a trajectory c(t) beginning at a 




(*») . 



dc 

~dt 



!F(c) • A K ■ Y . 



(1) 



positive vector c(0) = c° remains in the invariant polyhedron P :— (c° + S')n]R>Q. 
Multistationarity refers to the existence of more than one steady state in some 
invariant polyhedron. A chemical reaction network may admit multistationarity 
for all, some, or no choices of positive parameters fty. 

Horn initiated the investigation of small chemical reaction networks by enu- 
merating networks comprised of "short complexes," those whose corresponding 
monomials c v have degree at most two |14ll5j . Networks that consist of at most 
three short complexes do not permit multiple steady states. 

The next section establishes that the following graph, which we call the 
Square, depicts a smallest reversible multistationary chemical reaction network: 



«12 




«23 



«34 



In the horizontal reactions, two units of species one are transformed into two of 
species two (or vice- versa), while a third unit remains unchanged by the reaction. 
In the vertical reactions, only one is transformed. 

The Square appeared in non-reversible form as networks 7-3 in [16] and 4.2 
in [11]. The matrices whose product defines the dynamical system (fTJ) follow: 



*(<0 



CiC 2 , 



2- 



\ 



-K12 — K14, 

K21 




K12 

-K21 — K.23 







K23 
K43 



«14 



K34 

-K41 — K 43/ 



Y 



(Z 0\ 
1 2 
3 
\2 I J 

There may be two or even three steady states in each invariant polyhedron P; 
Example [1] in the next section provides a choice of positive rate constants Kij 
that give rise to three steady states. Sections 4 and 5 determine precisely which 
parameters give rise to two steady states and which yield three. Moreover, we 
compute this semi-algebraic parametrization for all networks on the same four 
vertices as the Square, in other words, networks with complexes cf, Cxc|, c|, 
and c\ C2- The parametrization is captured in Table 1 and can be computed "by 
hand," but larger systems may require techniques of computational real algebraic 
geometry pQ. For example, our problem of classifying parameters according to 
number of steady states is labeled as Problem P2 in [21], where it is addressed 
with computer algebra methods. 



3 The Smallest Multistationary Network 



Following equation (7) of [5], the Matrix- Tree Theorem defines the following 
polynomials in the rate constants of the Square: 



Theorem 7 of [5] provides an ideal Mg that is toric in these Ki coordinates, and 
the variety of Mq is the moduli space of toric dynamical systems on the Square. 
In this case, the ideal Mq is the twisted cubic curve in the Ki coordinates, 
generated by the 2x2-minors of the following matrix: 



Theorem 7 of [5] says that for a given choice of positive rate constants Kij, 
the equations |lj define a toric dynamical system if and only if the minors of 
the matrix ((21) vanish. In general the codimension of Mg is the deficiency of a 
network; see Theorem 9 of [5]. Here the deficiency is two. Recall that a toric 
dynamical system is a dynamical system ([l} for which the algebraic equations 
tf'(c) ■ A K — admit a strictly positive solution c* S R> > this solution is called a 
complex balancing steady state [16) . In this case there is a unique steady state in 
each invariant polyhedron P, so multistationarity is ruled out. Toric dynamical 
systems exhibit further nice properties; for details, see [5113116] . 

It is no coincidence that the original monomials of the Square, namely cf, 
C1C2, c 2 , c\c2, parametrize the twisted cubic curve. In fact, the following general 
result follows from Theorem 9 in [5] . 

Proposition 1. Assume that a chemical reaction network G is strongly con- 
nected and all monomials c Vi have the same total degree. Then the toric variety 
parametrized by tf'(c) coincides with the variety of Mg- 

For the Square, each one-dimensional invariant polyhedron P is defined by some 
positive concentration total T = c\ + c 2 . The steady states in P correspond 
precisely to the positive roots of the following cubic polynomial: 

ps(x) = (-2ki2 - k u )x 3 + (k 4 i - K 43 )a; 2 + (2k 2 i - k 23 )x + {k 32 + 2k 34 ) ; 

this polynomial arises by substituting x :— c\jc 2 into the equation dc\/dt = 
—dc2/dt. From this point of view, we reach some immediate conclusions. First, 
the algebraic degree of this system is three, which bounds the number of steady 
states. Second, the number of steady states and their stability depend only on 
the rate parameters Ky , and not on the invariant polyhedron P or equivalently 
the choice of total concentration T. Also, by noting that ps(x) is positive at 



Ki 
K 2 
K 3 



K23K34K41 + K21K34K41 + K21K32K41 + K21K32K43 i 

K14K32K43 + K 12«34K4i + K12K32K4! + K12K32K43 , 

K14K23K43 + K14K21K43 + K12K23K4I + K12K23K43 1 

K14K23K34 + K14K21K34 + K14K21K32 + K12K23K34 ■ 




K 2 
K 3 




(2) 



x = and is negative for large x, we see that the Square admits at least one 
steady state for any choice of rate constants. Recall that the discriminant of a 
univariate polynomial / is a polynomial that vanishes precisely when / has a 
multiple root over the complex numbers |20j . Maple computes the discriminant 
of ps to be the following polynomial: 

- 108k 12 K32 — 432k^ 2 ^32K34 — 432k^ 2 k 34 — IO8K12K14K32 

- 432K12K14K32K34 - 432K12K14K34 + 64^12^21 - 96«i 2 K2lK23 + 48^12^21^23 

- 72K12K21K32K41 + 144K12K21K32K43 - 144K12K21K34K41 + 288K12K21K34K43 

- 8K12K23 + 36K12K23K32K41 - 72K12K23K32K43 + 72K12K23K34K41 

- I44K12K23K34K43 - 27k 14 k 32 - IO8K14K32K34 - 108ki 4 k 34 + 32ki 4 k 21 

- 48Ki4K2iK23 + 24^14^21^23 _ 36ki4K21«32K41 + 72ki4K 2 iK32K43 

- 72K14K21K34K41 + 144K14K21K34K43 - 4K14K23 + I8K14K23K32K41 

- 36K14K23K32K43 + 36K14K23K34K41 - 72K14K23K34K43 + 4k| 1 k 41 

- 16K2lK4lK43 + 16^21^43 ~ 4k2iK23«4i + I6K21K23K41K43 - I6K21K23K43 
+ K 23 K 41 ~ 4K23K41K43 + 4k2 3 K 43 - 4k 32 k|i + 24k 3 2K 4 iK 4 3 - 48 k 3 2K41K 43 

+ 32k 32 k 43 - 8k 34 k|i + 48k 34 k 4 iK43 - 96K34K41K43 + 64k 3 4K 43 . 

As ps is cubic and has at least one positive root, its discriminant is negative 
if and only if p$ has one real root and one pair of complex conjugate roots; in 
this case, the Square has precisely one steady state. When the discriminant is 
non-negative, the system may admit one, two, or three steady states; we analyze 
this case fully in the next section. 

Example 1. Consider the following rate constants for the Square: 

(ki 2) K14, «2i, «23, K32, K34, «4i, K 43 ) = (1/4, 1/2, 1, 13, 1, 2, 8, 1) . 

This yields ps{x) = — x 3 + 6x 2 — 11 + 6, which has three positive roots: x — 1, 2, 
and 3. This is an instance of bistability; it is easy to determine that x = 1 and 
x = 3 correspond to stable steady states, while the third is unstable. In the next 
section we determine the conditions for an arbitrary vector of rate constants to 
admit one, two, or three steady states. 

Recalling the definitions given earlier, the Square has the following properties: 
the number of complexes is n = 4, the number of connected components of G 
is I = 1, the number of species is s — 2, and the dimension of any invariant 
polyhedron is a = 1. The main result of this section states that this network is 
minimal with respect to each of these four parameters. 

Theorem 1. The Square is a smallest multistationary, mass-preserving, re- 
versible chemical reaction network with respect to each of the following param- 
eters: the number of complexes, the number of connected components of G, the 
number of species, and the dimension of an invariant polyhedron. 



Proof. First I = 1 and a — 1 are clearly minimal. Next any mass-preserving 
system with n < 2 or s = 1 has no reactions or has deficiency zero. Finally, an 
n = 3 system has deficiency zero or one; in the deficiency one case, the Deficiency 
One Theorem of Feinberg rules out the possibility of multistationarity [T^] • □ 

Among all mass-preserving multistationary systems that share these four 
minimal parameters, the Square is distinguished because its monomials arc of 
minimal degree. A connected network of lower degree would consist of at most 
three of Horn's "short" complexes [2] . 

We now discuss the possible connection of the Square to biology by comparing 
it to the following simple network: 

C X Cy Cy (3) 

Cx * Cy . 

Network ([3} is a modified version of the following molecular switch mechanism 
proposed by Lisman fTS] : 

C-xCy > ('xy * Cy 

CyCp > Cyp > C x Cp. 

Here x denotes a kinase in an inactive state, y is the active version, and p is a 
phosphatase. In the first reactions, y catalyzes the phosphorylation of x, turning 
x into y; the second reactions correspond to dephosphorylation. By skipping 
the binding steps, making all reactions reversible, and noting that removing p 
effectively scales the second reaction rate constant, we obtain the network ©. 
The reactions of © are similar to c\c y c% and c\ ^» c§, which are reactions 
in the generalization of the Square network examined in the next section; this 
suggests the possible biological relevance of the reactions of the Square. For 
example c\ci — ► c\ can be viewed as a reaction in which species two catalyzes 
the reaction c\ — > c|. Such a positive feedback loop-in which a high amount of 
some species y encourages the further production of the same species-occurs in 
biological settings. For example, the recent work of Dentin et al. finds that high 
glucose levels in diabetic mice promote further glucose production in the liver, 
which is triggered by the binding of glucose (which we may view as y) to the 
transcription factor CREB (x) [5]. 

This paper focuses on the Square and more generally, the networks that share 
the same complexes as the Square. In the following section, we shall determine 
which of these are bistable. The one with the fewest edges is the only one with 
two connected components rather than one, and is featured in the last section. 

4 Parametrizing Multistationarity 

The aim of this section is similar to that of Conradi et al. [3] , which determined 
the full set of parameters that give rise to multistationarity for a biochemical 
model describing a single layer of a MAPK cascade. However we additionally 



determine the precise number of steady states: zero, one, two, or three, and 
determine their stability. The family of networks we consider are those that have 
the same four complexes as the Square. In other words, we classify subnetworks 
of the complete network depicted here: 




Each of the twelve rate constants By is permitted to be zero, which defines the 
parameter space R> 2 of dynamical systems. The main result of this section is 
summarized in Table 1, which is the semi- algebraic decomposition of the twelve- 
dimensional parameter space according to the number of steady states of the 
dynamical system. The conditions listed there make use of certain polynomials 
in the rate constants, including the signed coefficients of the polynomial p: 

50 = 2ni2 + 3ki 3 + K14 , 

51 = K41 — K42 — 2«4 3 , 

52 = — 2k 2 1 + K 2 3 — K 2 4 , 

53 = 3k 3 i + K 32 + 2k 34 , 
where p generalizes the polynomial ps from the Square: 

p(x) = -S a x 3 + Six 2 - S 2 x + S 3 . 



We now derive the entries of Table 1 for those networks without boundary 
steady states (this includes the case of the Square). These cases are precisely 
the ones in which So > and S3 > 0. Our approach is simply to determine the 
conditions on the coefficients of p for the polynomial to have one, two, or three 
positive roots. 

In this twelve-parameter case, the discriminant of p is a homogeneous degree- 
four polynomial with 113 terms. For the same reason as that for the Square, there 
is one steady state when the discriminant is negative. Now assume that the 
discriminant is non-negative. Then p has three real roots, counting multiplicity; 
recall that the positive ones correspond to the steady states of the chemical 
reaction network. Now the constant term of a monic cubic polynomial is the 
negative of the product of its roots, so by examining the sign of the constant 
term of p, we conclude that p has either one positive root and two negative roots, 
or three positive roots. Continuing the sign analysis with the other coefficients 
of p, we conclude that there are three positive roots if and only if S\ > and 
S*2 > 0. We proceed by distinguishing between the cases when the discriminant 
is positive or zero. If the discriminant is positive, then we have derived criteria 



Table 1. Classification of dynamical systems arising from non-trivial (having at 
least one reaction) networks with complexes c\, C\c\, c?,, c\ci- Listed are the 
number of steady states and the number of steady states that are stable. The 
discriminant of p is denoted by D. The signed coefficients of p are denoted by 
Sq, Si, S 2 , and 53. The triple root condition consists of the equations (4). 



Condition 


Steady states 


Stable states 


D < and So S3 = 








D < and else 


1 


1 


D > and S ,Si,S 2 ,S 3 > 


3 


2 


D > and S , Si, S 2 > and S :i = 


2 


1 


D > and Si, S 2 , S 3 > and S = 


2 


1 


D > and So = S3 = and S1S2 < 








D > and eZse 


1 


1 


D = and So, Si, S2, S3 > and triple root condition 


1 


1 


D = and So, Si, S2, S3 > without triple root condition 


2 


1 


D = and Si < So = < S 2 and S3 > 








D = and Si < S3 = < S 2 and So > 








D = and else 


2 


1 



for having one or three steady states; this is because the roots of p are distinct. 
If the discriminant is zero, then in the case of one positive root, the two negative 
roots come together (one steady state) . In the case of discriminant zero and three 
positive roots, then at least two roots come together (at most two steady states); 
a triple root occurs if and only if the following triple root condition holds: 

3S S 2 = Si and 27Sg£ 3 = Sf . (4) 

These equations are precisely what must hold in order for p to have the form 
p{x) = — (.x — a) 3 . Finally, stability analysis in this one-dimensional system is 
easy, and this completes the analysis for the networks without boundary steady 
states. The remaining cases can be classified similarly to complete the entries 
of Table 1. To parametrize the behavior of the Square, we simply reduce to the 
case when each of its parameters k 12 , K14, K21, K23, K 32, K 34, K41, and K43 are 
positive and all others are zero. 

By determining which sign vectors in (0, +) 12 can be realized by a vector 
of parameters that yields multistationarity, we find a necessary and sufficient 
condition for a directed graph on the four complexes of the Square to admit 
multistationarity. This condition is that the graph must include the edges labeled 
by K23 and K41 and at least one edge directed from the vertex cf or c\. In this 
case, for appropriate rate parameters arising from Table 1, the dynamical system 
has multiple steady states. Therefore, we can enumerate the reversible networks 
on the four complexes that admit multistationarity: there is one network with 
all six (bi-directional) edges, four with five edges, six (including the Square) 
with four edges, four with three edges, and one with two edges. These sixteen 



networks comprise the family of "smallest" multistationary networks. For the 
two-edge network, the decomposition from Table 1 is depicted in Figure [T] in the 
next section. 



5 Subnetworks of the Square 

Subnetworks of the Square are obtained by removing edges. From the parametriza- 
tion in the previous section, we know that up to symmetry between c\ and C2, 
only two reversible subnetworks of the Square exhibit multiple steady states. 

The first network is obtained by removing the bottom edge of the Square. In 
other words A K is replaced by 



A K 





V K 41 



«12 
-K2i - K 23 
K32 





«14 ^ 

K23 



-KuJ 



In this subnetwork, the four parameters of Theorem Q] are the same as those 
of the Square. The system is a toric dynamical system if and only if the following 
four binomial generators of Mq vanish: 



K14K32 
K12K32M1 
K14K21 



^23^41 , 

K14K21K23 
K12K41 , 



K12K32 — ^21^23 . 

We note that both K23 times the third binomial and K14 times the fourth bi- 
nomial are in the ideal generated by the first two binomials. Therefore, an 
assignment of positive parameters for this network defines a toric dynamical 
system if and only if the following two equations hold: K14K32 = K23M1 and 

K12K32K41 = K14K21K23- 

The second subnetwork of the Square is obtained by removing one additional 
edge, the one between the vertices labeled by c\ and Cic|. The new A K is 

f-K U K U \ 

-K 23 K23 
K32 -K32 

y K41 —kaiJ 

The network graph G is now disconnected, and p reduces to 

p{x) = - K14X 3 + k^x 2 — K23X + k 32 . 
The discriminant of p is 



D 



-27/C? 4 K32 ~ 4K14K23 + 18ki4K23«32K41 + «23 K 41 - 4/C32K4I 



Further, the toric condition reduces to the single equation 



^23^41 = K14K32 , 



which defines the Segre variety. A single equation suffices to define the space 
of toric dynamical systems; this corresponds to the fact that this subnetwork 
has deficiency one, while the previous subnetwork has deficiency two. The semi- 
algebraic decomposition of the previous section for this four-parameter network 
can be depicted in three dimensions by setting one parameter to be one, in other 
words, by scaling the equations ([I]); this is displayed in Figure [TJ 




Fig. 1. This depicts the semi-algebraic decomposition of Section 4 for the sub- 
network of the Square in which only the vertical edges remain and K41 — 1. At 
the left is the discriminant-zero locus. Parameter vectors lying below this sur- 
face give rise to dynamical systems with three steady states. Those above the 
surface yield one steady state; these include parameters of the toric dynamical 
systems, which are the points on the Segre variety which appears on the right. 
Parameters on the discriminant-zero locus correspond to systems with either one 
(if 3K14K32 = K23) or two steady states. This figure was created using Maple. 



We remark that Horn and Jackson performed the same parametrization for 
the following special rate constants: 

(«12, «14, «2ii «23, K32, k 34 , K41, k 43 ) = (e, 0, 1, 0, e, 0, 1, 0) , 

where e > 0. Their results are summarized as Table 1 in |16j . Their analysis 
notes that any instance of three steady states can be lifted to establish the same 
in the (reversible) Square. In other words, in a small neighborhood in H^>q of 
a vector of parameters that yields three steady states of the directed Square, 
there is a vector of parameters for the bi-directional Square that also exhibits 
multistationarity. The specific criterion for when lifting of this form is possible 
appears in Theorem 2 of Conradi et al. [I]. As this approach is widely applicable, 
further analysis of small networks may be fruitful for illuminating the dynamics 
of larger biochemical networks. 

We have seen that the family of Square networks is the smallest class of 
bistable mass-action kinetics networks. Whether nature has implemented one of 
these (perhaps with additional components to provide robustness) in a biological 
setting is as yet unknown, but it is also remarkable that these networks exhibit 
a simple switch mechanism, which we now explain. Consider the case of three 
steady states. The corresponding positive roots x\ < X2 < X3 of p in Section 4 
are the equilibria for the ratio of concentrations c\/c2. To switch from the low 
stable equilibrium x\ to the high stable equilibrium X3 is easy: simply increase 
the concentration ratio C\jc-i past X2, and the dynamics will do the rest. 
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